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ABSTRACT 

We apply a new method to determine the local disc matter and dark halo matter 
density to kinematic and position data for ~ 2000 K dwarf stars taken from the 
literature. Our method assumes only that the disc is locally in dynamical equilibrium, 
and that the 'tilt' term in the Jeans equations is small up to ~ 1 kpc above the plane. 
We present a new calculation of the photometric distances to the K dwarf stars, and 
use a Monte Carlo Markov Chain to marginalise over uncertainties in both the baryonic 
mass distribution, and the velocity and distance errors for each individual star. We 
perform a series of tests to demonstrate that our results are insensitive to plausible 
systematic errors in our distance calibration, and we show that our method recovers the 
correct answer from a dynamically evolved N-body simulation of the Milky Way. We 
find a local dark matter density of p dm = 0.025±g:oi3 M pc -3 (0.95±o;|g GeVcm" 3 ) at 
90% confidence assuming no correction for the non-flatness of the local rotation curve, 
and p dm = 0.022;° °^ M Q pc~ 3 (0.85±^ GeVcm" 3 ) if the correction is included. 
Our 90% lower bound on p dm is larger than the canonical value typically assumed in 
the literature, and is at mild tension with extrapolations from the rotation curve that 
assume a spherical halo. Our result can be explained by a larger normalisation for 
the local Milky Way rotation curve, an oblate dark matter halo, a local disc of dark 
matter, or some combination of these. 
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1 INTRODUCTION 

The local dark matter density is an average over a small 
volume, typically a few hundred parsecs, around the Sun. It 
provides constraints on the local halo shape and allows us 
to predict the flux of dark matter particles in laboratory de- 
tectors. The latter is required to extract information about 
the nature of a dark matter particle from such experiments, 
at least in the limit of a few tens to hundreds of detections 



( Peter 2011 1. The Galactic halo shape can be constrained by 



combining two methods of determining the local dark matter 
density. Firstly, one can infer it from the Galactic rotation 
curve (plm)- This requires an assumption about the shape of 
the Galactic halo (typi cally spherical; e. g.|Sofue et al.|2009| 



Weber fe de Boer|2010||Catena fc Ullio|2010[ ). Secondly, one 
can calculate the dark matter density locally from the ver- 
tical kinematics of stars near the Sun (pdm) (e. g. |Bahcall| 
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1984b| |Holmberg fe Flynn||2000| ). If p dm < pf£, this sug- 
gests a prolate dark matter halo for the Milky Way; while 
Pdm > Pdmi could imply either an oblate halo or a dark disc 



( Lake|1989||Read et al.|2008| [Readet al.|2009 i. 

Determining the local matter density from the kinemat- 
ics of stars in the Solar Neighbourhood has a long history 
dating back to[Oort]([l932|[l^0| in the 1930's. Oort used the 
classical method of solving the combined Poisson-Boltzmann 
equations for a sample of stars, assumed to be stationary in 
the total matter distribution of the disc. He found 50% more 
mass than the sum of known components. A more modern 
study by Bahcall (1984b) introduced a new method that 



described the visible matter as a sum of isothermal compo- 
nents. He also found dyn amically significant dark matter in 
the dis(p"| ( [Bahcall 1984a). Using faint K dwarfs at the South 



1 We should be careful about what wc mean by 'dark matter 
in the disc'. Early studies likejOort (1932J were typically inter- 
ested in missing disc-like matter (a 'thin dark disc'); more mod- 
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Galactic Pole, Bahcall et al. ( 1992 1 confirmed his earlier re- 



sult that more than 50% of the mass was dark, although with 
a lower statistical significance. However, the early studies 
by [Oortl fl932"l [l960| and |Bahcall| ( |1984b|a| ) assumed that 
different tracers could be simply averaged to form a sin- 



gle tracer population. Kuijken & Gilmore (1989a) demon- 



strated that the two samples of F stars analysed by [Bahcall] 
( 1984b I were not compatible (i.e. they had different spatial 



density distribution, but no evidence for a difference in their 
kinematics) and therefore should not be averaged. They re- 
analysed the K giant sample used by Bahcall (1984a I, as- 



signing more realistic errors to the density profile and using 
a more detailed fit to the velocity data, finding a value of to- 
tal matter density compatible with the observed one. They 
concluded that the determination of the local volume den- 
sity remained limited by systematic and random errors with 
the available data. 

With the launch of the ESA satellite Hipparcos (1997), 
the kinematics and position of tracer stars were measured 
with much higher accuracy. The improved distance measures 
give a much more accurate measurement of the local lumi- 
nosity function, so that the total amount of visible mat- 
ter can be better estimated as well. The latest dynamical 
measurements of the local density of matter - p to t - from 
Hipparcos data show no compelling evidence for a signifi- 
cant amount of dark matter in the disc ( Creze et al.||199"8 



Holmberg fc Flynn||2000[ ) . |Holmberg fc Flynn| ( |2000| | found 
Ptot = 0.102 ± O.OIMqpc - ' 3 , with a contribution of about 
O.O95M0pc~ 3 in visible matter, consistent with the Kuijken 



fc Gilmore| |1989|'s value. 

In addition to the local volume density, several authors 
have calculated the local surface density of gravitating mat- 
ter, probing up to larger heights above the disc plane (typ- 



ically ~ 1 kpc; e.g. Kuijken & Gilmore 


1989b| |Kuijken &| 


Gilmore 1989 1991| |Holmberg & Flynn 


20041. Using faint 


K dwarfs at the South Galactic Pole, and using a prior from 


the rotation curve, Kuijken & Gilmore 


(1989 19911 find 



from the rotation curve assuming a spherical Galactic dark 
matter halep] (e. g. Sofue et al 



2009; Weber & dc Boer 2010 



Catena fc Ullio|2010 1. A similar result was found in the post- 



Hipparcos era by Holmberg & Flynn ( 2004 1 . Recently, Moni 



Bidin et al. (20121 have estimated the surface density using 
tracers at heights 1.5 < z < 4 kpc above the disc, making 
a rather stronger claim (incompatible with the earlier re- 



sults of Kuijken & Gilmore ( 1991 1 and Holmberg & Flynn 



(2004 1 ) that there is no dark matter near the Sun. However, 
Bovy & Tremaine (2012 I demonstrate that this result is er- 



roneous and owes to one of ten assumptions used by |Moni] 



Bidin et al. ( 2012J) being false. Furthermore, |Sanders| ( |2012[ ) 
estimate that the velocity dispersion gradients derived by 



ern studies try to constrain a significantly more extended dark 
matter halo that has a near-constant dark matter density up to 
~ 1 kpc. Even the 'dark disc' predicted by recent cosmological 
simulations ( Read et al. 2008 Read et al. 2009 ) is sufficiently hot 
that its dark matter distribution is approximately constant up to 
~ 1 kpc. Throughout this paper when we talk about 'dark matter 
in the disc' we refer to a constant density dark matter component 
within the disc volume. 

2 Note that this consistency with the rotation curve is somewhat 
circular since this is input as a prior in their analysis. 



Moni Bidin et al. ( 2012 1 could be biased by up to a factor of 
two, which would also significantly alter their determination 

of Pdm- 

With next generation surveys round the corner (e.g. 
Gaia; Jordan 2008), a significant improvement in the num- 
ber of precision astrometric, photometric and spectroscopic 
measurements is expected. For this reason, |Garbari et al.| 
(20111 (hereafter Paper I) revisited the systematic errors in 
determining pd m from Solar Neighbourhood stars; these will 
likely soon become the dominant source of error, if they are 
not already. We were the first to use a high resolution N- 
body simulation of an isolated Milky Way-like galaxy to gen- 
erate mock data. We used these mock data to study a pop- 
ular class of mass modelling methods in the literature that 
fit an assumed distribution function to a set of stellar trac- 
ers ( |Holmberg fc Flynn|200"o[|Binney fc Tremaine|2008l ) . We 
found that realistic mixing of stars due to the formation of a 
bar and spiral arms (similar to those observed in the Milky 
Way) breaks the usual assumption that the distribution 
function is separable, leading to systematic bias in the recov- 
ery of pdm- We then introduced a new method that avoids 
this assumption by fitting instead moments of the distribu- 
tion function (i.e. that solves the Jeans-Poisson equations). 
Our Minimal Assumption method (or MA method) uses a 
Monte Carlo Markov Chain technique (hereafter MCMC) to 
marginalise over remaining model and measurement uncer- 
tainties. Given sufficiently good data, we showed that our 
method recovers the correct local dark matter density even 
in the face of disc inhomogeneities, non-isothermal tracers 
and a non-separable distribution function. 

In this article, we apply our MA method to real data 
from the literature. The key advantages of our new method 
over previous works are that: (i) we use a 'minimal' set of 
assumptions; (ii) we use a MCMC to marginalise over both 
model and measurement uncertainties; and (iii) we require 
no prior from the Milky Way rotation curve as has been 
commonly used in previous works (Kuijken fc Gilmore|1989| 



1991 Holmberg fc Flynn 2000 2004 1. This latter means that 



we can compare our determination to that derived from the 
rotation curve to constrain the Milky Way halo shape. Our 
method requires at least one equilibrium stellar tracer pop- 
ulation with known density fall off v{z) and vertical velocity 
dispersion a\{£), both as a function of height z. The require- 
ments for a suitable sample of stellar tracers are that: (i) 
they are in dynamical equilibrium with the Galactic poten- 
tial (i.e. they must be sufficiently dynamically old to have 
completed many vertical oscillations through the Galactic 
plane); (ii) they are available in sufficient numbers to give 
good statistical precision; (iii) they have reliable distances 
and vertical velocities v z ; (iv) the sample completeness needs 
to be sufficiently well understood in order to measure the 
density fall off as a function of the distance z from the disc 
plane; and (v) they extend up to 2-3 times the disc scale 
height (in order to break a degeneracy between the disc and 
dark matter densities; see Paper I). While full six dimen- 
sional phase space information is now available for a large 
number of stars (e.g. RAVE Steinmetz 2003 , Stcinmetz et al. 
[2006] [Zwitter et al.|[2008] and SEGUE |Yanny et al.||2009[ ), 
these surveys are magnitude rather than volume complete, 
with additional survey selection effects based on colour. This 
makes it difficult to reliably estimate v(z) for a given tracer 
population. For this reason, we return to the volume com- 
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plete K-dwarf data from Kuijken & Gilmore ( 1989 1 for our 
disc tracers - the 'KG' data. These data consist of a photo- 
metric sample of 2016 K dwarf stars, complete in the z-range 
~ 0.2 — 1.5 kpc, with a spectroscopic sample of 580 K dwarfs 
(most of which are included in the photometric catalogue). 
We use data from Hipparcos and SEGUE (Kotoneva et al. 
|2002"} |Zhang||2012"l ) to perform a new photometric distance 
measurement for each K dwarf star. We model the local 
gravitational potential using the baryonic mass distribution 
of the Galactic disc by |Flynn et al. (20061. 

This article is organised as follows. In Section 2] wc 
present the K dwarf data from | Kuijken fc Gilmore ( 1989 1 
(hereafter KG989II) and describe our new distance deter- 
minations 1 2.1 1. In Section [3] we summarise our MA mass 
modelling method (3.1 1 and, for comparison, the method 
adopted by KG89II (3.2 1. In Section 3.3 we test both meth- 



ods on a mock data set derived from a dynamically evolved 
N-body simulation. In Section [4] we apply our MA method 
to the KG data and present our results. Finally in Section 
[5] we summarise and present our conclusions. 
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2 DATA 

KG89II present a catalogue of 2016 K stars with photome- 
try in B and V bands, and another of 580 K dwarfs (most 
of which are also included in the photometric catalogue) in- 
cluding radial velocities, at South Galactic pole (figure [T]). 
The stellar density fall off of these tracers was derived from 
star counts. At large z, the mean metallicity of the stars 
is known to decrease below the Solar Neighbourhood value. 
Such a gradient translates into an absolute magnitude gra- 
dient, since the position of the main sequence in the colour- 
magnitude diagram changes with the metallicity: metal poor 
dwarfs are fainter than metal rich ones at the same temper- 
ature or colour (the opposite is true for giant stars). So, if 
there is a vertical metallicity gradient, the photometric par- 
allaxes used for the derivation of the density fall off - v{z) 
- will be systematically wrong as one moves away from the 
plane. Unfortunately, the metallicity could not be measured 
directly for these stars, so KG89II derived the density fall 
off using an assumed constant metallicity gradient for the K 
dwarfs. They considered two different gradients to estimate 
the magnitude of the uncertainties, namely <f[Fe/H]/dz = 
(constant metallicity) and d[Fe/~H]/dz — — 0.3dexkpc _1 (in 
their analysis KG89II consider this latter as the fiducial 
metallicity gradient for K dwarfs). 

The tracers' vertical distance determination is funda- 
mental for our analysis. Twenty years after KG89II's study, 
we can re-calibrate the distances for these stars using mod- 
ern survey data to estimate the metallicity distribution func- 
tion of K dwarfs at different z, and Hipparcos parallaxes 
to calibrate the photometric distances. Our distance re- 
calibration procedure is described, next. 



2.1 A new distance determination for the K 
dwarfs 

To calculate the vertical distances z for KG89II's sample, we 
must derive a relationship between the metallicity [Fe/H], 
the vertical distance z and the absolute magnitude My of K 
dwarf stars in the disc. Since the metallicity is not included 



Figure 1. Spatial distribution of the K dwarf sample from 
KG89II. Blue crosses: spectroscopic sample; empty circles: pho- 
tometric sample. The red dashed lines mark the range of 200 < 
z < 1200 pc use in our analysis (see Section [2TTJ . 



in KG89II's catalogue, we can only hope to derive a distance 
distribution function P, (z) for each star of the sample, based 
on an observed metallicity distribution function for these 
stars. 

We consider two different catalogues of K dwarfs with 
distances and metallicity for our calibration. The first cata- 



logue by Kotoneva et al. (2002) consists of 431 K dwarfs from 



Hipparcos, representing a complete catalogue of the metal 
content in nearby K dwarfs extending up to z ~ 100 pc. The 
vertical distances z for these stars are very accurately de- 
termined from Hipparcos parallaxes. The second catalogue 



by |Zhang| ( |2012| > contains 5000 SEGUE K dwarfs spanning 
a much wider range of z, namely between 300 and 2000 pc. 
However, in this case, the distance determination for these 
stars is much less certain: the distance errors are about 10%. 

We combine these two catalogues to build the metal- 
licity distribution function (MDF) Q([Fe/H](z), z) for K 
dwarfs. Comparing the SEGUE metallicity distribution 
function at z = 500 pc with the MDF from Kotoneva et al.| 
( 2002 ) for z < 100 pc (see the red and black solid histograms 
in figure [5] respectively), we notice that the two MDFs have 
very similar shape, but shifted means. There is a vertical 
metallicity gradient from to 500 pc of ~ — 0.4dexkpc _1 
(corresponding to a shift of — 0.2dex; see figure [5] black dot- 
ted histogram). At larger height than this, the gradient is 
weaker: the SEGUE MDF at 1 kpc (dashed red histogram) 
is similar to the one at 500 pc. We adopt the Kotonev a et al.| 
( 2002 1 MDF in the Galactic plane, then we apply a linear 
metallicity gradient of ~ — 0.4dexkpc _1 between 100 pc and 
500 pc to match the SEGUE MDF at z = 500 pc; we extend 



this shifted |Kotoneva et al.| ( |2002[ )'s MDF up to z = 800 pc, 
and we adopt the SEGUE MDF for z > 800 pc, as shown in 
figure [3] We explore an alternative Q([Fe/H](«), z) for the 
K dwarfs in Appendix [A] 

With the metallicity distribution function 
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Figure 2. The z < 100 pc MDF (black solid histogram) from |Ko-| 
|toneva et al.| (2002fl . A shift of 0.2 dex (black dotted histogram) 
approximately overlaps the MDF from SEGUE K dwarfs com- 
puted at z ~ 500 pc (red solid histogram). The red dashed his- 
togram is the SEGUE MDF at z ~ 1000 pc. 
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Figure 3. Co mbined MDF of |Kotoneva et al. | j2002| ) and SEGUE 
(Zhang 2012]) . The colours show the probability values of [Fe/H] 
given z. 



Q([Fe/H](z), z) for K dwarfs, we next derive the dis- 
tance probability distribution P*(z) for each star of the 
KG89II's catalogue. This requires calculating the relation- 
ship between ^-distance and metallicity [Fe/H] for each 
star. 

Figure [4] show s the absolute magnitude My (ordinate) 
of Kotoneva et al. (2002 1 's K dwarfs as a function of colour 



index B — V (abscissa) and metallicity [Fe/H] (colours). We 
fit this using a polynomial: 

M v =a ,o + ai,o(S - V") + a ,i[Fe/H] + a 2 ,o(B - V) 2 

+ 01,1 (B - V0[Fe/H] + a , 2 [Fe/H] 2 + a , 3 (B - V") 3 
+ 03,i (S - V) 2 [Fe/R] + ai, 2 (B - V)[Fe/R] 2 

(1) 

The best-fit parameters are [-5.795, 27.92, 0.1291, -22.74, 
-2.003, 0.04917, 7.113, 1.02, -0.04274], with an error of 
~ 0.03 mag. 

Once we have M v = M V (B - V, [Fe/H]), we write the 
distance modulus as: 

V-A v -M v +5 

d = 10 s = d(V, B-V, [Fe/H]) (2) 

where V is the apparent magnitude and Ay is the extinction; 
we use the value Av = 0.062 mag from Schlcgcl et al.| ( |T998[ ), 
given the mean Galactic coordinates of KG89II's data. The 
vertical distance z for a single star is then: 



z,(l,b,d) = z*(l, b,B — V, [Fe/H]) 



(3) 



where I and b are the Galactic longitude and latitude. We 
know I, b and B~V for each star of KG89II's sample, so the 
only free parameter is [Fe/H]. This means that the vertical 
distance for each star will be given by a probability dis- 
tribution P*(z) corresponding to a metallicity distribution 
P«([Fe/H]) for that star: 



P^) = P*([Fe/H](z)) 



(4) 



In practice, this equation must be solved iteratively because 
[Fe/H] is itself a function of z through equation [3] The iter- 
ative process proceeds as follows: 

(i) We start the first iteration by assuming that the dis- 
tance distribution function of a single star is P*(z) = 1 for 
all the possible z([Fe/H]) calculated using equation[3] 

(ii) The MDF marginalised over z for a single star is given 
by: 



P,([Fe/H])= / Q([Fe/E](z),z)P,(z)dz 
Jo 



(5) 



where <2([Fe/H](z), z) is the observed MDF at each z for all 
K dwarfs as described previously. 

(iii) A new P»(z) is calculated through equation [4] using 
the P*([Fe/H]) just computed. 

(iv) We restart steps (ii) 



to 



(iii) to calculate P*([Fe/H]) 



with the new P*(z) until it converges. For all stars, we ob- 
tained convergence in less than 5 iterations. 

The density fall off of the photometric sample and the 
velocity dispersion of the spectroscopic one, obtained with 
our new distance estimates, are shown in figure [5] For our 
analysis, we use the density and the velocity dispersion pro- 
files only over the range 200 < z < 1200 pc (red dashed 
lines). This assures that our sample is volume complete and 
avoids significant contamination by K giant stars Kuijkcn & 
Gilmore| ( |1989[ ). The corresponding quantities computed by 
KG89II, assuming a metallicity gradient of — 0.3dex/kpc, 
are plotted as a comparison (empty circles). Our new den- 
sity profile and velocity dispersion do not differ greatly from 
those of KG89II; they are compatible within the quoted er- 
rors. In Appendix [XJ we explore the effect of a rather ex- 
treme variation in the assumed MDF - ignoring |Kotoneva] 
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et al. ( 2002 I and SEGUE data - finding that our results are 



not sensitive to plausible changes in our distance calibra- 
tion. The referee of this article pointed out that the study 
of Kotoneva et al. ( 2002[ ) has been updated by Casagrande 
|et al.| ( |2007[ ). The two studies are very much compatible, 
but the scatter in equation [T] of 0.03 mag becomes 0.27 mag 
when the newer data are used. We tested the impact of this 
larger scatter in magnitude on the distance calibration, find- 
ing that the density and velocity dispersion profile remain 
unchanged in the range of z of interest, with a negligible 
increase in the uncertainties (see Appendix [A]). 



3 METHOD 

3.1 The MA method 

The MA method presented in Paper I uses the Poisson- Jeans 
system to predict the density fall off of a tracer population 
in a given gravitational potential. The comparison between 
this predicted density fall off and the observed one allows us 
to constrain the gravitational potential and, consequently, 
the underlying dark matter distribution. 

Here we summarise the basic equations; for a detailed 
description of the MA method see Section 2.1 of Paper I. 

The MA method is based on three main assumptions: 

(i) The system is in equilibrium (steady state assump- 
tion). 

(ii) The dark matter density is constant over the range of 
\z\ considered. 

(iii) The 'tilt' term (iiVcrfj z ) in the cylindrical Jeans 
equation: 

- M{ Rua Rz ) + -^a z )+u l -=0 (6) 
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Figure 5. Upper panel: K dwarf stellar density profile (filled cir- 
cles with error bars) derived from a Monte Carlo sampling of 
P*(z) for each star. As a comparison, the density profile (assum- 
ing a metallicity gradient of —0.3 dex/kpc) from KG89II is plotted 
as empty circles. Lower panel: The similarly derived vertical ve- 
locity dispersion as a function of z (filled circles with error bars). 
The corresponding determination by KG89II is represented by the 
empty circles. In both panels, the two red dashed lines show the 
range of z considered in our analysis, over which the photometric 
sample is volume complete and we avoid significant contamina- 
tion from K giants stars. 



is negligible compared to all other terms. Here v, a z and a\ z 
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are the number density and the velocity dispersion compo- 
nents of a tracer population moving in potential $. 

With these assumptions, the Jeans equation becomes a func- 
tion only of 2 and we can neglect the other two Jeans equa- 
tions in R and 9: 



oz 



d® do* 
1 

dz dz 



0. 



(7) 



Solving this equation for a single tracer population, we ob- 
tain its density v(z) at each height z: 



Kg) 
v{zo) 



Q 2 (z ) 



exp 



1 d$ 



al(z) dz 



dz 



(8) 



Given the density at the midplane p s ,j(0) and the vertical 
velocity dispersion o 2 j(z) as a function of z for each of the 
gas and stellar populations in the local disc, we can model 
the full disc density distribution as a superposition of such 
elements: 



a i,j dz 



dz 



(9) 



In Paper I, we showed that accurate measurement of the ver- 
tical velocity dispersion of the tracers o 2 (z) is crucial, how- 
ever in the mass modelling we can assume that all the visible 
matter components are isothermal - i.e. o 2 j = o\ j (0) - and 
equation [9] simplifies to: 



(10) 



The Poisson equation then determines the potential $ from 
the density. In cylindrical polar coordinates, this is given by: 



d 2 ® 

dz 2 



AnG(p s (z) + pt 



with: 



Pdm = Pdm(R) - {A-kGRY 



OR 



(11) 



(12) 



where pdm(-R) is the halo mass density (following assumption 
|(ii)| this is assumed to be independent of z in the volume 
considered); and V C (R) = (Rd$>/dR) 1/2 is the (total) circu- 
lar velocity at a distance R (in the plane) from the centre of 
the Galaxy. For a flat rotation curve, the second term van- 
ishes and Pdm(-R) = Pdm(R)- The rotation curve correction 



can be calculated from the Oort constants .4 and B (Binney 
|fc Merrifield|1998[ ): 



(4itGR) 



-idV 2 



B 2 -A 2 



(13) 



OR 2-kG 

We solve equations [10] and [II] numerically for a given 
tracer population, adopting the following procedure: 

(i) We make initial trial guesses for p B ,j(0), pdm, and the 
vertical velocity dispersion for the visible matter component 
in the plane o 2 j(0). 



(ii) We solve equation 10 to obtain &(z) and its first 



derivative §§, with $(0) =~|f| = 0. 

(iii) We insert this result into equation [8] for the vertical 
density fall off of the tracers v p (z) and we compare this with 
the observed distribution v(z) to obtain a goodness of fit. 



This procedure requires many input parameters, such as the 
normalisations and dispersions of each of the disc compo- 
nents, and the vertical dispersion profile of the tracers (typ- 
ically poorly constrained). To explore the parameter space 
and marginalise over the uncertainties, we use a Monte Carlo 
Markov Chain (MCMC) method. 

In Paper I we showed that our method is able to recover 
the correct value of the local dark matter density, even in 
presence of large visible matter density fluctuations due to 
the spiral arms. Notice that the MA requires no prior from 
the Milky Way rotation curve, as has been commonly used 
in previous works; this means that we can compare our de- 
termination to that derived from the rotation curve to con- 
strain the Milky Way halo shape. 



3.1.1 Application of the MA method to real data 

When we apply the MA method to real data, we must deal 
with distance and velocity uncertainties, and account for 
survey geometry and/or the sample completeness. In par- 
ticular, the K dwarf data from KG89II are not assigned dis- 
tances, but instead, as described in Section [2.1[ a z proba- 
bility distribution function P, (z) for each star of the sample. 
In addition, the vertical velocities for each star are measured 
with an uncertainty of about 0.5 — lkm/s. To marginalise 
over these uncertainties, we proceed in the following way: 

(i) For each model n of the MCMC, we select a different 
vertical distance z* for each star in the sample, according to 
its z distribution function P*(z). We make sure that for each 
star included in both the spectroscopic and the photomet- 
ric sample we pick a unique z^ value. For each star of the 
spectroscopic sample we also draw a vertical velocity value 
v* n from a Gaussian distribution, according to its velocity 
error bar. 

(ii) We bin the data in z to construct the observed density 
fall off v n {z) and velocity dispersion o 2 n (z) for the tracers, 
selecting stars with z between 0.2 and 1.2 kpc. The velocity 
dispersion is calculated using the velocity scale algorithm 



described in Beers et al. ( 1990 ) 



(iii) We use the trial guesses for p SJ ,n(0), o z ^ n (0) and 
Pdm,n to solve equation |10| and |11| to obtain $> n (z) an d its 
first derivative ^$ IL . 

OZ 

(iv) We insert this result and the vertical velocity disper- 
sion of the tracers ol„(z) into equation [8] to predict the 
tracers' density fall off. When applying the MA method to 
the real data, we add the visible matter surface density E s 
as a further constraint. For each model of the MCMC we 
compute S s as 



Es.Tl 2 



/ ps(z)dz = I fh,j,n{°) ex P 2 

JO JO a V CT z,j,n 



(14) 

(v) We compare the predicted density fall off v n ,p{z) with 
v n (z), the predicted visible matter surface density Sf n with 
the observed one (E s ) and we calculate the corresponding 
Xn, accepting or rejecting the model n. 

(vi) When a model is accepted, we restart from |(i)| with 
the following model (n + 1) and so on, exploring the whole 
parameter space. 
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3.2 The KG method 

KG89II used the K dwarf data to calculate the total surface 
density at the Sun position E(-Rq). Their approach (the KG 



method) is similar to the method adopted later by Holmberg 
fc Flynn|p000||2004[) (the HF method), that we analysed in 



detail in Paper I. 

Instead of measuring the vertical velocity dispersion of 
the tracers as a function of z to predict the density fall 
off, the HF method uses the tracers' velocity distribution 
function in the mid-plane of the Galactic disc / z (v z ,o); as- 
suming that the distribution function is separable, / = 
f R ,e(v R ,v e ,R) x f z (v z ,z). [HoTmberg fc Fly£n] pOO) [2TJ04 1 
integrate this distribution function over z-velocities to pre- 
dict the density fall off of the tracers (see Section 2.2 of 
Paper I for more details): 



/oo p 
f z (v z , z)dv z = 2 / 
-OO J $ 



f z (E z )dE z 



(15) 



where f z (z,v z ) = f z (E z ) and E : 
energy. This equation can be written as: 

fz{Vzfi)v z fidv z fl 

- 2$(z) 



* W ^/2[E Z - $(z)] 

\v\ + &(z) is the vertical 



(16) 



where v z o is the vertical velocity at the Galactic mid-plane 
(z = 0).' 

The MA demands only that the tilt term in the Jeans 
equation |6| is small with respect to the other terms, the HF 
method requires the stronger assumption that the z-motion 
(and so the distribution function) is completely separable 
from the motion in the radial and azimuthal directions; this 
latter implies that the tilt term is exactly zero. 

The HF approach has the advantage of exploiting the 
whole available information about the shape of the veloc- 
ity distribution function of the tracers. However, we demon- 
strated in Paper I that when the separability of the distribu- 
tion function is not fulfilled, the HF method leads to biased 
results. Using a high resolution simulation of a Milky Way 
like galaxy, we showed that the onset of spiral arms and a 
bar can cause significant radial mixing that breaks the sep- 
arability of the motion in the z and R directions, violating 
this key assumption of the HF method. This effect becomes 
increasingly important with height z. It is important to no- 
tice that the HF method can not be corrected for this bias 
since the separability of the potential (and of the distribu- 
tion function) lies at the heart of the method. By contrast, 
in the MA method, the separability of the potential enters 
only in the neglected tilt term of the Jeans equation (that is 
assumed to be small as compared to the other terms). If the 
tilt term is for some reason large - i.e. the radial derivative 
of the density weighted tilt of the velocity ellipsoid is large 
- a correction can straightforwardly be applied to our MA 
method. However, we expect the tilt term to be small (see 
Paper I and Binney & Tremaine 20081, and so we do not 



consider such a correction in this paper. 

KG89II's approach relies on the same key assumption 
about the distribution function as the HF method. In most 
studies, the density v(z) of the tracers is known to better 
precision than the velocity distribution f z (v z ,z). For this 
reason, KG89II work in the opposite direction with respect 
to |Holmberg fc Flynn|p000l[2C]0ll ) and predict f z (v z ,z) from 
the observed v(z). Applying an inverse Abel transform to 



equation |15| they obtain: 
1 



-dv/d$ 



Kj E ; - E z 



(17) 



so there is a unique relation between v(<f>) and f z (E z ). No- 
tice that f z (E z ) depends on v(&(z)) only at large z, where 
the potential exceeds E z , i.e. beyond z = $ _1 (_E Z ). Thus 
an additional key advantage of the KG method is that one 
can model the potential at large distances from the Galactic 
plane, ignoring the detailed distribution of matter at small 
z. KG89II parameterised the gravitational potential $(2) 
above the bulk of the disc matter (where it is sensitive only 
to the total surface density of gravitating matter) as: 



$(«) = K(^z 2 + D' 2 - D) + Fz 2 



(18) 



where D is the disc scale height, K is proportional to the 
total disc surface density E(_Rq), and F oc (the effective 
halo density). KG89II used a range of Galactic mass models 
(calculated using different values of the disc mass M, the 
radial disc scale-length R^, the circular velocity V c (Rq) and 
Sun position Rq) to ensure consistency with the Galactic 
rotation curve (assuming a spherical Milky Way halo) and 
therefore to obtain a relation between F and K. Note that 
already this is different from our MA approach where we 
use no information about the rotation curve to constrain 
our mass models. 

Given the observed space density of a tracer population 
v(z) and a set of gravitational potential models <&(z), one 
can solve equation |17| To reduce the noise in the differential 
of v(z), KG89II fitted it with a double exponential. KG89II 
then used the derived f z (E z ) for each potential model $(z) 
to compute the likelihood of the spectroscopic sample: 



a 



f °° f z (E z )dE z 



(19) 



where the product is over all stars in the spectroscopic sam- 
ple, and select the potential parameters that maximise this 
likelihood function C. 

The KG method, like the HF method, uses the full 
shape of the observed velocity distribution function, max- 
imising the use of the available information. It is also con- 
venient because it does not require a detailed model of 
the gravitational potential or an accurately measured tracer 
density fall off close to the Galactic plane. However, its draw- 
back is that, like the HF method, it relies on a key assump- 
tion that the vertical distribution function is only a function 
of E z . In the following section we test, using the high res- 
olution N-body simulation described in Paper I, how this 
assumption affects the result derived using the KG method. 

3.3 Testing the methods using an N-body 
simulation 

Before applying the MA method to the real K dwarf data, 
we use the most dynamically evolved stage of the simulation 
described in Paper I as a mock data set, to test the effect 
of the velocity errors and of the asymmetric distance distri- 
bution function P*(z) on the MA method's result. Then we 
use the same mock data to probe how the non-separability 
of the tracers' distribution function affects the KG method. 
The mock data consist of a high resolution (30 million 
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Figure 6. Results for recovering p^m an d p s from the simulation. Top left panel: The recovered dark matter density (red filled circles) 
for all 8 volumes analysed around the disc (90% and 68% confidence intervals are shown). The actual value of the dark matter density 
is marked as a blue filled circle. Top right panel: The density of models explored by the MCMC projected onto Pdm _ Ps space for the 90° 
volume. The blue dot shows the true value for p^ m and p s ; the diagonal cyan blue line shows the total matter density; and the red dot 
shows the median recovered p<j m and p s . with 90% errors marked. Bottom left panel: Histogram of the recovered P(j m from the MCMC 
mode ensemble for this same volume. The striped grey area is the 90% confidence interval. The cyan dashed line and shaded area give 
the actual value of P(j m with error bars. Bottom right panel: Testing the effect of assuming a separable distribution function. Each dot 
in this plot shows the likelihood L n of each MCMC model calculated assuming a separable distribution function (see equation |17| the 
plot shows — log(£ n ), so the more likely models have a lower ordinate value). The red dashed line shows a fit to the points. Notice that 
assuming that the distribution function is separable produces a bias towards low p<j m (the true distribution function for the simulation 
is not separable; see figure u\. This bias effect was observed in all the volumes considered. 



disc star particles and 15 million halo dark matter particles) 
N-body simulation of an isolated Milky Way like galaxy. 
The initial conditions were built to contain some thousand 
stars in the volume size required for our analysis. With the 
dynamical evolution of the simulation, the disc developed a 
bar and spiral arms. For more details about the simulation 
features and how it compares to the real Milky Way, see 
Section 3.1 of Paper I. 

We consider several volumes in the disc of the simu- 
lated galaxy at a distance of Rq = 8.5 kpc from the cen- 
tre. The MA method solves the Jeans equation for a one- 
dimentional slab (equation [7p, so the radial size of the vol- 
umes, AR = 0.25 kpc, is chosen to fulfil this approximation, 



but still contain an enough large number of stars (see Section 
3.3.1.1 of Paper I for more details). 

As described in Section [3. 1.1| we assign a different ve- 
locity value v* n and a different z* to each star at every iter- 
ation n of the MCMC. In the application of the MA method 
to the simulation, the velocity values are drawn from a Gaus- 
sian distribution centred on the true velocity value and with 
a width of lkm/s, while the z values are selected from a 
lognormal distribution around the true value. Because of 
the numerical resolution of the simulation, we cannot fit 
the density profile up to 1.2 kpc, since at such height we 
quickly run out of star particles and the velocity dispersion 
is poorly measured, so we use stars with 0.2 ^ z ^ 0.75 kpc. 
For the simulation data, we model the visible mass distribu- 
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tion as a single component, characterised by its mass den- 
sity p s ,j = ps (0) and its velocity dispersion on the mid- 



plane 



erj(0). We let the dark matter density freely 



vary between and 0.2 M© pc , and the other parame- 
ters - p s (0) and a z (0) - vary inside their error bars. We 
adopted Poisson errors for the velocity dispersion <r 2 (0) and 
the current uncertainties on the total visible matt er density 
in the plane (i.e. 0.014 M Q pc" 3 



4.1[ ) for p s (0). 
which can 



see Section 

We include the rotation curve correction term, 
be easily computed for the simulation, in the calculation. 
We test convergence of our MCMC chains by starting with 
pdm seeded at two different values (namely pdm = and 
Pdm = O.2M0 pc~ 3 ) and running until the two chains are 
statistically indistinguishable. 

In figure [5] the results for the MA method are shown. 
The upper left panel shows the results for eight different vol- 
umes around the simulated disc. Notice that in all cases, the 
mean correct answer (blue filled circles) is recovered within 
the 90% confidence interval, while for four out of eight of the 
patches (with a fifth at 45° extremely close) it is recovered 
with the 68% confidence interval. This is consistent with our 
confidence intervals having the meaning of a purely statis- 
tical error, despite each patch being systematically different 
(each patch samples a different region of the disc with differ- 
ent local dynamics). The remaining three panels focus on the 
results for the volume at 90° . The top right panel shows the 
density of models explored by the MCMC projected onto 
Pdm-Ps space; the bottom left panel shows a histogram of 
the dark matter density for all the models explored by the 
MCMC; and the bottom right panel explores the effect of 
assuming a separable distribution function (of which more, 
next). 

In Paper I, we noticed that, at this evolved stage of 
the simulation, the distribution function f(vR,ve,v z ,R,z) 
of the tracers is not separable. This leads to the HF method 
producing biased results - either underestimating or over- 
estimating the local dark matter density (see Section 3 of 
Paper I for more details). 

The separability of the distribution function lies at 
the heart of the KG method too. To reproduce a KG- 
like method, we consider all v n {z) and model potentials 
$„(«) explored by our MA method MCMC chain. For each 
model in the chain, we then calculate a distribution function 
fz,n{E z (v z ,&)) through equation 17 and use it to compute 
the likelihood C n of the velocity data via equation |19| 

In practice, the integral in the denominator of equa- 
tion [l9] is calculated numerically from Ef in = 1 to _B" mx = 
7000 km 2 s -2 , which is chosen to avoid the divergence at 
E z — and to ensure we cover all energies of interest (the 
contribution of the high energy tail of f z ,n(E z ) is negligible 
with respect to the low energy part; see figure [7b. 

In the bottom right panel of figure [6j the likelihoods 
C n of the velocity data are plotted against the correspond- 
ing values of the local dark matter density for the differ- 
ent MCMC models pd m ,n- This panel shows that there is an 
anti-correlation between the computed likelihood £.„ and the 
corresponding value of the local dark matter density pdm,n'- 
the likelihood of the velocity data is larger (i.e. — log(£ n ) 
is lower) for gravitational potential models with low pdm,™; 
this means that we expect the KG method to artificially 
favour low dark matter density values. For all the explored 
volumes around the disc, we always obtain this same anti- 
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Figure 7. The distribution function — / - of the tracers in the 
simulation. The black line shows distribution function measured 
directly from the simulation averaged over a volume 8.375 < R < 
8.625 kpc; 0.2 < z < 0.75 kpc. The orange dashed line shows 
fz(E z ) predicted from the density fall off v(z) via equation |17| 
assuming the correct value of p E and pdm in the computation of 
the gravitational potential. The blue and red dashed lines show 
fz(Ez) predicted from u(z), assuming pdm = and twice Pdm, 
respectively. All distribution functions are normalised by the in- 
tegral of fz(Ez) between the minimum and the maximum E z of 
the star particles in the volume considered. 



correlation, so the bias on pd m will have always the same 
sign. 

To understand this effect, in figure [7] we show the dis- 
tribution function f z (E z ) calculated from the same density 
profile v{z), but using three potentials with different values 
of Pdm, namely pdm = (blue dashed line), the true value of 
Pdm (orange dashed line) and twice the true pdm (red dashed 
line); the black line represents the actual distribution func- 
tion of the stars in the volume. The potential corresponding 
to the true value of the dark matter density does not predict 
the distribution function correctly, while with pdm = we 
obtain a better agreement with the measured f z (E z ). It is 
clear from this plot that the likelihood of models with low 
dark matter density, calculated through equation |19[ will be 
higher than the true model. Therefore the KG method will 
be biased towards low pam- 



4 RESULTS 

4.1 Measuring the local matter and dark matter 
density 

In the previous section, we showed that the MA method 
is able to recover the correct dark matter density within 
our 90% confidence interval, even in presence of asymmetric 
distance errors, velocity uncertainties and a non-separable 
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distribution function. We now apply the MA method to 
the real K dwarf data, proceeding as described in Section 
|3~1~T1 The mass distribution in the Galactic disc is mod- 
elled as a superposition of 15 isothermal components, listed 
in Table [I] As parameters to fit in the MCMC, we use 
the local dark matter density pdm, the total visible den- 
sity in the midplane p s (0), and the relative fractions of 
the visible components p St j(0) / p s (0) and their velocity dis- 
persions in the midplane a z ,j- We allow the densities and 
the velocity dispersions of the different components to vary 
within their measured uncertainties (the errors for each 
component are given in Table [T]). We let the total visible 
density in the plane p s (0) vary within its observed range: 
p s (0) = 0.0914 ± 0.0140 M Q pc" 3 (extrapolated from Table 
[T|; and we let the dark matter density vary between and 
O.2M pc~ 3 . For each model explored by the MCMC, we 
calculated the visible surface density Ef jn through equation 
[14] and compare it with the total surface density we obtain 
from Table [l] i.e. E° bs ± AE° bs = 49.4 ± 4.6 M Q pc" 2 , to 
calculate the total \ 2 '- 



2 _ 2 
X Xsurf 



where: 



and 



Xsurf 



(S° 



2 

Xu 



yp 



(AE° 



X, 



= E 



(20) 



(21) 



(22) 



where the sum is extended to all the bins and L\u n i are the 
uncertainties on the density fall off. 

The rotation curve correction term can be calculated 
from the Oort constants A and B. To determine the Oort 
constants, we must use stellar tracers that are well-mixed. 
The most recent estimates of A and B from F giants 



( Branham| |2010[ ) and K-M giants (|Mignardj I2000J) from 
Hipparcos give A = 14.85 ± 7.47 km s~ kpc - and B — 
-lO.SSie.SSkms" 1 kpc" 1 and A = 14. 5± 1.0 km s" 1 kpc" 1 
and B — —11.5 ± 1.0 kms -1 kpc -1 , respectively. Averaging 
these two values we obtain a correction term of —0.0033 ± 
0.0050 M Q pc" 3 . We test for convergence of the MCMC by 
starting several chains at different initial values of all the 
parameters and running until they are statistically equiva- 
lent (after removing an initial burn-in phase of 100 accepted 
models for each chain). 

The results for the MA method applied to the real 
K dwarf data are shown in figure [8] The upper panel 
shows the density of the models explored by the MCMC 
(grey contours) in the p s — pdm plane; the median with 
90% errors is shown by the black dot and corresponds 



to p dm = 0.025±g;gi| M pc" 



(0.95Z 



IGeVcnT 



'0 



adding the rotation curve correction, we obtain pdm = 



M pc" 



(0.85 



+0.57 



Ge V cm , see the red dot in 



figure |8| . 

We expect the MA method to primarily constrain the 
total matter density in the plane p E + pdm, which means 
that we should see an oblique degeneracy between p B and 
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Figure 8. Upper Panel: The recovered visible and dark matter 
densities. The grey contours are the density of models explored 
by the MCMC. The black dot shows the median recovered value 
of p s and pdm with 90% errors marked; the red dot shows the 
same but including a correction for the local non-flatness of the 
Milky Way rotation curve (-0.0033±0.0050 M pc -3 ). The pur- 
ple area represents the values estimated by Holmberg Sz Flynn 
(12000 p (this appears as a diagonal stripe since these authors only 
constrained pd m + p s ) . The blue-dashed lines show our priors on 
p a . The horizontal dashed orange line marks the Standard Halo 
Model value of p dm (0.008 M pc -3 ). The green area marks the 
range of p^m we extrapolate from KG91. Lower Panel: The re- 
covered dark matter density as a function of the visible matter 
surface density S s . The meaning of the symbols is the same as in 
the upper panel. The striped grey area is the range in the visible 
matter surface density determined by |Flynn et al.| l ]2006ft . 
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Table 1. The disc mass model taken from Flynn et al. 2006. For 
each component in the table, we give the local mass density in 
the midplane p(0) in M pc _ 3 ; the total column density £ in 
M pc — 2 ; and the vertical velocity dispersion <r z j(0) in kms -1 . 
Uncertainties on the densities are assumed to be 50% for all the 
gas components (indicated with *) and 20% for all of the stellar 
components. The largest uncertainties come from the gas that 
remains poorly constrained (compare, for example, compilations 
in Flynn et al. 2006, ? and |Ferriere| j2001| ). For the thick disc, 
the column density is well known, while the velocity dispersion 
and the volume density are poorly known such that they should 
have larger error bars. However, these two quantities are essen- 
tially nuisance parameters for our analysis here. Since they anti- 
correlate and - as pointed out by |Kuijken &: Gilmore| ( |1989b[ ) — 
the local gravitational potential is mainly constrained by the col- 
umn density, we simply assume small errors for both here such 
that the integrated column agrees with the observed value. 



Component 


fi,o(0) 




^(O) 1 / 2 




[Ms pc" 3 ] 


[M pc- 2 ] 


[kms" 1 ] 


H* 


0.021 


3.0 


4.0 ± 1.0 


HI(1)* 


0.016 


4.1 


7.0 ± 1.0 


HI(2)* 


0.012 


4.1 


9.0 ± 1.0 


Warm gas* 


0.0009 


2.0 


40.0 ±1.0 


Giants 


0.0006 


0.4 


20.0 ±2.0 


My < 2.5 


0.0031 


0.9 


7.5 ±2.0 


2.5 < M v < 3.0 


0.0015 


0.6 


10.5 ± 2.0 


3.0 < M v < 4.0 


0.0020 


1.1 


14.0 ±2.0 


4.0 < M v < 5.0 


0.0022 


1.7 


18.0 ±2.0 


5.0 < My < 8.0 


0.007 


5.7 


18.5 ± 2.0 


My > 8.0 


0.0135 


10.9 


18.5 ± 2.0 


White dwarfs 


0.006 


5.4 


20.0 ±5.0 


Brown dwarfs 


0.002 


1.8 


20.0 ±5.0 


Thick disc 


0.0035 


7.0 


37.0 ±5.0 


Stellar halo 


0.0001 


0.6 


100.0 ± 10.0 



Pdm in the p B — pd m plane, as was observed for the simula- 
tion (see figure Rjb. However, unlike the simulation that has 
only one visible matter component in the disc, our real-data 
mass model comprises some 15 separate components with 
different scale heights. This introduces new freedoms that 
wash out the oblique degeneracy in the p s — pdm plane. If 
we plot instead, however, the total surface density of visible 
matter in the disc (bottom panel of figure |8| the degener- 
acy is once again clearly visible as a diagonal elongation of 
the MCMC model density contours. Notice that many of 
the models explored by the MCMC lie outside of the range 
given by the Flynn et al. 2006 mass model (grey striped 
band). This simply means that the prior we placed on E a is 
not very strong. However, the full area explored is in good 
agreement with the more conservative measurements of the 
total visible matter surface density by Kuijkcn & Gilmorc 
( 199l|) (hereafter KG91), namely E s = 48 ±8 M pc" 2 ; and 
bylFlynn & Fuchsl dl994|), E s = 49 ± 9 M pc" 2 . Even if we 



include only models that lie within the grey striped region, 
our recovered pdm is little affected. 

In Appendix |Bj we test the robustness of our result, 
exploring how it changes if one considers either only the 
low z bins (0.2 < z < 0.7 kpc) or only the high z bins 
(0.6 < z < 1.2 kpc). We find that the low z data do not 
provide any information about pdm, but they still favour 
low E s compared to the prior we imposed. The low z bins 
present a noisier and not monotonically increasing velocity 



dispersion, however they are better constrained and domi- 
nate the x 2 nt. Using the high z data, we lose information 
about the disc and E s settles into the centre of its prior dis- 
tribution; this leads to a systematically lower pdm, since the 
sum of the two is well constrained. This tells us that the 
origin of our high pd m is the slightly lower E s required by 
the velocity dispersion data near the plane. 

We also tested the effect of changing the assumed errors 
on the baryonic mass model. We first reduce them to an op- 
timistic 10% error for the stellar normalisations and 30% for 
the gas normalisations; this has little effect on the resulting 
determination of pd m because sufficient freedom remains in 
our mass model to allow degeneracies between pd m and the 
baryonic components (this reflects the broken degeneracy 
seen in figure |8|. We then increased the errors by removing 
the prior on p s (0) altogether. This also has a small effect on 
our recovered pd m - This confirms the results shown in Ap- 
pendix[B]that it is the velocity dispersion data, not our prior 
that are constraining our mass model. The low z (< 500 pc) 
data constrain the surface density profile of the disc, while 
the high z data (> 500 pc) are dynamically sensitive to dark 
matter. 

The green horizontal stripe in the upper panel of figure 
[8] marks the value of pdm we extrapolate from KG89II and 
KG91. KG91, using the same K dwarf data analysed in this 
paper, determine the total dynamical surface density up to 
1.1 kpc: E dyn = 71 ±6 M pc- 2 . If we subtract from this 
the contribution of the observed visible matter E s = 48 ± 
8 M pc~ 2 , we can calculate pdm, assumed to be constant 
in the range < z < 1.2 kpc, as: 



KG 
Pdm 



Jdyn 



2 • 1100 



(23) 



This gives: = 0.010 ± 0.005 M pc" 3 . 

Our new result from our MA method is in tension with 
pJm, obtained from the same data set. This could owe ei- 
ther to our different distance calibration or to our new MA 
method that does not require any assumption about the 



separability of the distribution function (see Section 3.31. In 
figure [5j we already showed that our new distance calibra- 
tion does not significantly affect the velocity dispersion and 
the density fall off of the K dwarfs. In addition, in Appendix 
|A"} we explore the effect of using a constant metallicity gradi- 
ent for the K-dwarfs of — 0.3dexkpc" , exactly as assumed 
by KG89II. This metallicity distribution is not compatible 
with modern data; we only use it to illustrate the sensitiv- 
ity of our results to the MDF of the K dwarf stars, and 
to fully understand why our determination of pdm is larger 
than that of KG91. Using the KG89II's MDF, our recovered 
value of pd m is slightly smaller and therefore in better agree- 
ment with KG89II and KG91. However, our median recov- 
ered value remains significantly larger than the upper bound 
of the KG91 result. This suggests that our new distance de- 
terminations are not the primary reason for the systematic 
shift. 



In our tests on the N-body simulation (Section 3.3 1, we 
showed that, when the distribution function of the tracers 
is not separable, the method adopted by KG89II leads to a 
systematic underestimate of pa m . In the lower panel of fig- 
ure^ we plot the likelihood C„ of each model explored by 
our MCMC - calculated through equation [19] - against the 
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corresponding value of pd m - Unlike the similar plot for our 
simulation data (figure[6] bottom panels), there is now a sig- 
nificant vertical dispersion in the models. This owes to the 
increased freedom present in our 15-parameter mass model 
for the real-data. However, the highest likelihood models 
(the bottom envelope of points in the plot) show a similar 
trend as seen for the simulation data: higher likelihood mod- 
els have systematically smaller pdm. We conclude that the 
primary difference for our larger value of pdm as compared 
to KG89II is that our MA method requires no assumption 
about the separability of the distribution function. 

In the upper panel of figure [9] we plot a histogram of 
Pdm from all the models explored by the MCMC. The striped 
area is the 90% confidence interval (corresponding to the 
black dot of figure [8|; the result including the rotation curve 
correction is shown by the red error bar. Notice that our 
90% lower bound is larger than the Standard Halo ModeQ 
typically assumed in the literature (marked by the vertical 
dashed orange line). For a comparison, we plot the ranges 
of pdm at the solar radius obtained by Iocco et al. (20111, 



combining microlensing and rotation curve measurements, 
and using different halo models: the blue error bar corre- 
sponds to a spherical halo, while the cyan and purple bars 
correspond to oblates halos with potential flattening q = 0.9 
and q = 0.7, respectively. The magenta bar represents the 
dark matter density in presence of a dark disc, contributing 
0.25 — 1.5 times the dark matter (spherical) halo density, as 
predicted by |Read et al.] ( |2009 1 . 

From figure |9j we can see that our recovered density 
is in mild tension with the result for a spherical Milky Way 
halo. Moving to an oblate halo significantly reduces this ten- 
sion, however a flattening of q — 0.7 is likely inconsistent 
with measurements of the halo shape from the Sagittarius 
stream of stars (e.g. Ibata et al.|2001 1. If we wish to explain 
our median value for pd m that is very much larger than the 
canonical SHM value assumed in the literature, we require 
a local disc of dark matter that raises pd m without signifi- 
cantly altering the rotation curve. Interestingly, our median 
value is in excellent agreement with the range of dark discs 



predicted for our Galaxy by Read et al. (20081; Read et al 
(20091. 
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5 DISCUSSION AND CONCLUSIONS 

We have presented a new measurement of the local matter 
and dark matter densities from the kinematics of K dwarf 
stars near the Sun. We presented a new photometric dis- 
tance calibration for the the K dwarf data of KG89II (the 
KG data), derived using modern survey catalogues and the 
Hipparcos satellite data. We then used these data as tracers 
of the local gravitational potential to calculate the visible 
(p s ) and dark matter (pdm) densities at the solar position 
Rq and the surface density of the Milky Way disc up to 
1.1 kpc above the plane (S s ). 

To determine pd m and p s , we applied our new mass 
modelling method (presented already in Paper I) that re- 
lies on a minimum set of assumptions (the MA method) to 

4 The SHM is an isothermal sphere model for the Milky Way's 
dark matter halo with a value of the dark matter velocity disper- 
sion assumed to be Cj so ~ 270kms _1 . 
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Figure 9. Upper panel: A histogram of the recovered pd m from 
our MCMC chains for the MA method applied to the real K dwarf 
data. The striped grey area is the 90% confidence interval. The 
orange dashed line is the SHM value of pd m ; the blue, cyan and 
purple error bars correspond to the value of Pd m obtained by |Iocco] 
|et al.| ( |20 1 1 [ ) from a combination of microlensing and rotation 
curve data, with a spherical halo (potential flattening q = 1) 
and two oblate halos, with q = 0.9 and q = 0.7, respectively. 
The magenta error bar is the value of pdm expected if the Milky 
Way has a dark disc contributing 0.25 — 1.5 times the density of 
the (spherical) halo. The red error bar corresponds to our result 
after adding the rotation curve correction. Lower panel: The effect 
of assuming a separable distribution function. Each dot shows 
the likelihood of a given MCMC model calculated assuming a 
separable distribution function. Notice that the assumption of 
separability biases the result towards low Pd m - The orange and 
the green dashed lines have the same meaning as in the upper 
panel of|8] 
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the rejuvenated KG data. The key advantages of our new 
method are that: (i) we do not require any hypothesis about 
the shape of the tracers' velocity distribution function; (ii) 
we use a MCMC to marginalise over uncertainties in the 
distances and velocities of the tracer stars, and the under- 
lying baryonic mass model for the visible disc; and (iii) we 
require no prior from the Milky Way rotation curve as has 
been commonly used in previous works. This latter means 
that we can compare our determination to that derived from 
the rotation curve to constrain the Milky Way halo shape. 
We used a dynamically evolved high resolution N-body sim- 
ulation of a Milky Way-like galaxy as a mock data set to 
test our MA method, finding that we could correctly re- 
cover pdm and p s within our 90% confidence interval (for 
eight sample Solar neighbourhood-like volumes) even in the 
face of disc inhomogeneities, non- isothermal tracers, asym- 
metric distance errors and a non-separable tracer distribu- 
tion function. Furthermore, we confirmed the result from 
our Paper I that assuming a separable distribution function 
(as has been typically done in the modern literature) leads 
to a biased determination of pdm- 

Applying our MA method to the K dwarf data, we ob- 
tain a new measurement of the local dark matter density: 
p dm = 0.025lg;^3 M pc" 3 (0.95l^;^GeVcm- 3 ) ; which, 
adding a correction for the local non-flatness of the ro- 
tation c urve correction term (~ —0.0033 ± 0.0050, see 
Sections 3.1 and 4k, gi 
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larger than the results from KG89II and KG91 derived from 
the same data. We show that this primarily owes to our new 
MA modelling method (and the fact that it does not assume 
a separable distribution function for the tracers); our new 
distance determination for the K dwarfs plays a more minor 
role. At the same time we determine a value of the local visi- 
ble matter density of p s = 0.098±J5;{5?1 M pc" 3 that largely 
reflects the prior from our baryonic mass model. 

Our error bars are larger than is often quoted in 
the literature, however they reflect the full combination of 
model systematic, measurement and statistical uncertain- 
ties. Other recent determinations either rely on the rotation 
curve and therefore a strong assumption about the Milky 
Way halo shape, or require a large number of assumptions 
with associated (and typically unmodelled) systematic er- 
rors. 

In addition to measuring pdm and p s , we also obtain 
an estimate of the baryonic disc mass up to z = 1.1 kpc 
above the disc plane: E s = 45.5!ts'9 Mqpc -2 at 90% con- 
fidence. This is slightly lower than the mean of the prior 
from our baryonic mass model: E v i s = 49.4 ± 4.6 Mq pc -2 . 
Splitting the number into the contribution from stars and 
stellar remnants: E* = 33. 4+52 M© pc -2 and gas: E g = 
12. OOI2I5 Mqpc -2 , we see that our model favours slightly 
lower surface density in both the gas and the stars than 
the mean of our priors (E° bs = 13.3 ± 3.4 M Q pc" 2 , EJ bs = 
36.1 ±3.0 Mq pc" 2 ). 

It is this tendency for our models to favour lower disc 
surface density that leads to our high median value for pdm 
(see figure [8] and Appendix [B| . Unfortunately, current es- 
timates of the stellar and gaseous inventory in the Solar 
neighbourhood are too uncertain to confirm or rule out our 
Bovy et al.|2012| |Ferriere|20"oT l. 



favoured E s (e.g. 



Our median value of the local dark matter density is 



larger at 90% confidence than the Standard Halo Model 
value of pl m M = 0.008 M Q pc -3 (0.30 GeVcm -3 ) usually 
adopted in the literature. This could be a statistical fluc- 
tuation; in one out of eight patches in our simulated mock 
data, our method overestimated pdm by ~ 90%. However, if 
our high median value is confirmed by future data then it 
has some interesting implications. Firstly, it is particularly 
important for direct detection experiments because it im- 
plies a larger flux of dark matter particles and therefore a 
greater chance of detection. Secondly, our result is at mild 
tension with the value of p d m extrapolated from the rotation 
curve measurements, assuming a spherical dark matter halo. 
This suggests that the halo of our Galaxy is oblate and/or 
that we have a disc of dark matter, as predicted by recent 
cosmological simulations (see upper panel of figure [9| . 
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APPENDIX A: TESTING THE ROBUSTNESS 
OF THE MDF 

The most uncertain quantity in our re-analysis of KG89II's 
data is the variation of the K dwarfs' metallicity distribu- 
tion function with z, Q([Fe/H](z), z). In this appendix, we 
investigate how the adopted metallicity distribution func- 
tion affects the result of our analysis by exploring a differ- 
ent model for this function. We use the gradient adopted by 
KG89II, i.e. — 0.3dexkpc _1 , and set the mean metallicity to 
at z = (see upper left panel of figure [Al|. This MDF 
is not compatible with modern metallicity data for the K 
dwarfs (see figure we use it simply to illustrate our sen- 
sitivity to the assumed MDF, and to aid comparisons with 
the earlier KG89II results. 

we show the velocity dispersion (upper 



In figure Al 



right panel) and the density fall off (lower left panel) of the 
ttracers derived using the above MDF. In the lower right 
panel of figure |A1| the recovered dark and visible matter 
density are shown. The value of pdm obtained is slightly 
lower than that derived using our default MDF, namely 
Pdm = 0.022lS;Sii M pc- 3 , or p dm = 0.018±g;gi| M pc" 3 
including the rotation curve correction. However, our me- 
dian value for is still high and the overall result remains 
in tension with the SHM value (even when using this in- 



correct MDF). We conclude that our results are robust to 
plausible variations in the assumed MDF. 

In addition to the above, we also tested the impact of a 
different extinction value used in equation [2] on our results. 
We chose Ay = and Ay = 0.1. For Ay = (Ay = 0.1) the 
distances are slightly overestimated (underestimated) with 
respect to the extinction values considered in this article. 
This does not affect much the density fall off, but it trans- 
lates mainly in a slightly flatter (steeper) velocity dispersion. 
This leads to a small decrease (increase) of the recovered 
value of pdm- However the impact of the extinction is very 
small compared to the previous test shown in this Appendix. 

Finally we tested the impact on our distance calibration 
of the larger scatter in equation [T] obtained using th e more 
modern K dwarf catalogue from Casagrande et al. (20071 
instead of Kotoneva et al. (2002 I. The two studies are very 



much compatible, but the scatter in equation [T] of 0.03 mag 
increases to 0.27 mag when the newer data are considered. 
We used the 104 K dwarf stars from Casagrande e t al.| ( [2007 \ 
to build the relationship between My and (B — V , [Fe/H]), 
similar to equation [l] but using the polynomial down to 
o 2 , i (B - V) 2 [Fe/H] . Unlike |Kotoneva et al.| ( |2002[ > , the new 
catalogue also provides the error on the parallax, therefore 
we use a x 2 instead of a simple linear- least-square fit. To 
account for the 0.27 mag uncertainties on My, we convolved 
the distance probability distribution function with a Gaus- 
sian of width a — 13% (corresponding to 0.27 mag in My). 
The resulting density and velocity dispersion profiles are un- 
changed in the range of z of interest, with only a negligible 
increase in the uncertainties. Unlike in the case of the ex- 
tinction test, where distances were (slightly) systematically 
overestimated or underestimated, the magnitude scatter has 
only the effect of slightly increasing the random errors. The 
effect is weak so as long as a sufficiently large number of stars 
is used per bin, and the scatter is not too large. For this rea- 
son our determination of the local dark matter density is 
not affected by the increased scatter found by Casagrande 



et al. (20071 
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Figure Al. Left upper panel: The MDF adopted by KG89II: a constant metallicity gradient with height of 0.3dexkpc —1 normalised 
to a metallicity of 0. dex at z = kpc. Right upper panel: Stellar density profile derived from Monte Carlo sampling the probability 
distribution of z (dots with error bars). As a comparison, the density profile from KG89II is plotted as empty circles. The two red dashed 
lines show the completeness range. Lower left panel: The vertical velocity dispersion as a function of z. The black dots with error bar 
are derived from the probability distribution of z. Lower right panel: The projection of MCMC models onto the pdm-Ps plane for our 
MA method applied to the K dwarf data, using the MDF shown in the top left panel to determine the distances. The colours and the 
symbols are as in figure [8] 



APPENDIX B: EXPLORING THE 
ROBUSTNESS OF OUR p dm DETERMINATION 

In this Appendix, we explore the robustness of our determi- 
nation of pdm by analysing a low z (0.2 < z < 0.7 kpc) and 
a high z (0.6 < z < 1.2 kpc) subset of the KG data. The 
results are shown in Figure |B1| If we consider only the low 
z data (left panel), there is no information about the dark 
matter density. However, the data still favour low E s com- 
pared to our prior. If only the high z data are used (right 
panel), we lose the information about the disc and E s set- 
tles into more or less the centre of its prior distribution. This 



leads to a systematically lower /9dm- The above suggests that 
the origin of our high median pdm is the lower E s favoured 
by the velocity dispersion data close to the plane. 
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Figure Bl. The recovered dark matter density as a function of the visible matter surface density S s . Left panel: considering only low z 
bins (0.2 < z < 0.7kpc); central panel: considering the full z range: 0.2 < z < 1.2kpc (as lower panel of figurc|s]l; right panel: considering 
only high z bins (0.6 < z < 1.2 kpc). The meaning of the symbols is the same as in figure [8] The striped grey area is the range in the 
visible matter surface density determined by |FIynn et al.| l |2006| , that we used as a (weak) prior. 



